#library(calculus)

library(haven)
library(MASS)
library(Deriv)
library(xtable)

##############################################
#Choose EWMD, OMD, DWMD
weighting_matrix=2     #1=EWMD, 2=DWMD, 3=OMD#
##############################################


R_dataset <- read_dta("R:/SharedProjects/Shared2020-070/2016/extend_to_2020/JPE_Replication_dta/R_dataset_mom_se_dwmd.dta")

m=R_dataset$M1
V=data.matrix(R_dataset)

m=V[1:13,1:1]
alpha1=V[14,1]
var_alpha1=V[14,15]

V=V[1:13,2:14]

m
V


W=solve(diag(diag(V)))

for(i in 1:length(m)) {
  nam=paste("m",i,sep="")
  assign(nam,m[i])
}

f=function(b) {
  sigma_omega=b[1];sigma_v=b[2];sigma_phi=b[3];zeta=b[4];mu=b[5];delta=b[6];theta=b[7];lamda=b[8];sigma_psi=b[9];psi=b[10];sigma_csi=b[11];sigma_hrs=0
  mt1=(zeta-mu)/sqrt(1+sigma_omega^2+sigma_psi^2)
  mt2=(zeta-mu-delta)/sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2+sigma_hrs^2)
  mt3=(1+sigma_omega^2+sigma_psi^2)/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2+sigma_hrs^2))
  mt4=-(mu-theta)/sqrt(1+sigma_omega^2+sigma_psi^2)
  mt5=(zeta)/sqrt(1+sigma_omega^2+sigma_psi^2)
  mt6=sigma_psi^2/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_omega^2+sigma_psi^2))
  mt7=-1/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_csi^2))
  mt8=lamda*zeta+psi
  mt9=sqrt(lamda^2+sigma_phi^2)
  mt10=lamda/sqrt(1+sigma_omega^2+sigma_psi^2)
  mt11=lamda/sqrt(1+sigma_csi^2)
  mt12=lamda/sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2+sigma_hrs^2)
  mt13=-1/(sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2+sigma_hrs^2)*sqrt(1+sigma_csi^2))
  m_vector=c(mt1-m1,mt2-m2,mt3-m3,mt4-m4,mt5-m5,mt6-m6,mt7-m7,mt8-m8,mt9-m9,mt10-m10,mt11-m11,mt12-m12,mt13-m13)
  t(m_vector)%*%W%*%m_vector
}

beta=nlm(f,c(1,1,1,1,1,1,1,1,1,1,1))
b=beta$estimate

for(i in 1:length(b)) {
  nam=paste("b",i,sep="")
  assign(nam,b[i])
}

sigma_omega=b[1]; sigma_v=b[2];sigma_phi=b[3];pi=b[4];gamma=b[5];tau=b[6]
theta=b[7];lamda=b[8];sigma_psi=b[9];psi=b[10];sigma_csi=b[11]

b_to_print=rbind(pi,gamma,tau,theta,lamda,psi,sigma_omega,sigma_csi,sigma_psi,sigma_v,sigma_phi)
b_to_print=rbind(b_to_print,alpha1)
b_to_print



  f1=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{(zeta-mu)/sqrt(1+sigma_omega^2+sigma_psi^2)}
  f2=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{(zeta-mu-delta)/sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2)}
  f3=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{(1+sigma_omega^2+sigma_psi^2)/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2))}
  f4=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{-(mu-theta)/sqrt(1+sigma_omega^2+sigma_psi^2)}
  f5=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{(zeta)/sqrt(1+sigma_omega^2+sigma_psi^2)}
  f6=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{sigma_psi^2/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_omega^2+sigma_psi^2))}
  f7=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{-1/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_csi^2))}
  f8=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{lamda*zeta+psi}
  f9=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{sqrt(lamda^2+sigma_phi^2)}
  f10=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{lamda/sqrt(1+sigma_omega^2+sigma_psi^2)}
  f11=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{lamda/sqrt(1+sigma_csi^2)}
  f12=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{lamda/sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2)}
  f13=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{-1/(sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2)*sqrt(1+sigma_csi^2))}

fmn1=f1(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn2=f2(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn3=f3(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn4=f4(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn5=f5(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn6=f6(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn7=f7(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn8=f8(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn9=f9(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn10=f10(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn11=f11(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn12=f12(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)
fmn13=f13(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11)



fmn=rbind(fmn1,fmn2,fmn3,fmn4,fmn5,fmn6,fmn7,fmn8,fmn9,fmn10,fmn11,fmn12,fmn13)
fmn



ff1=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi) 
{(zeta-mu)/sqrt(1+sigma_omega^2+sigma_psi^2)}
ff2=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{(zeta-mu-delta)/sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2)}
ff3=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{(1+sigma_omega^2+sigma_psi^2)/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2))}
ff4=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{-(mu-theta)/sqrt(1+sigma_omega^2+sigma_psi^2)}
ff5=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{(zeta)/sqrt(1+sigma_omega^2+sigma_psi^2)}
ff6=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{sigma_psi^2/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_omega^2+sigma_psi^2))}
ff7=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{-1/(sqrt(1+sigma_omega^2+sigma_psi^2)*sqrt(1+sigma_csi^2))}
ff8=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{lamda*zeta+psi}
ff9=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{sqrt(lamda^2+sigma_phi^2)}
ff10=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{lamda/sqrt(1+sigma_omega^2+sigma_psi^2)}
ff11=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{lamda/sqrt(1+sigma_csi^2)}
ff12=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{lamda/sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2)}
ff13=function(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
{-1/(sqrt(1+sigma_omega^2+sigma_psi^2+sigma_v^2)*sqrt(1+sigma_csi^2))}

g1=Deriv(ff1)
g2=Deriv(ff2)
g3=Deriv(ff3)
g4=Deriv(ff4)
g5=Deriv(ff5)
g6=Deriv(ff6)
g7=Deriv(ff7)
g8=Deriv(ff8)
g9=Deriv(ff9)
g10=Deriv(ff10)
g11=Deriv(ff11)
g12=Deriv(ff12)
g13=Deriv(ff13)

zeta=pi;mu=gamma;delta=tau

h1=g1(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h2=g2(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h3=g3(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h4=g4(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h5=g5(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h6=g6(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h7=g7(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h8=g8(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h9=g9(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h10=g10(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h11=g11(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h12=g12(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)
h13=g13(sigma_omega,sigma_v,sigma_phi,zeta,mu,delta,theta,lamda,sigma_psi,psi,sigma_csi)

dfdb=rbind(h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11,h12,h13)

z=t(dfdb) %*% W %*% dfdb
covth=solve(z) %*% t(dfdb) %*% W %*% V %*% W %*% dfdb %*% solve(z)

se=sqrt(diag(covth))

for(i in 1:length(se)) {
  nam=paste("se",i,sep="")
  assign(nam,se[i])
}

se_to_print=rbind(se4,se5,se6,se7,se8,se10,se1,se11,se9,se2,se3,sqrt(var_alpha1))
se_to_print
test_stat=b_to_print/se_to_print
pval=2*pnorm(-abs(test_stat))

results=cbind(round(b_to_print,digits=4),round(se_to_print,digits=4),round(abs(test_stat),digits=2),round(pval,digits=4))

diffr=m-fmn
A=dfdb%*%z%*%t(dfdb)
P=diag(nrow=length(diffr))-dfdb%*%solve(z)%*%t(dfdb)%*%W
C=ginv(P%*%V%*%t(P))

mf=t(diffr)%*%C%*%diffr
dof=length(diffr)-length(b)
pvaluechi=pchisq(mf,dof,ncp=0,lower.tail=FALSE,log.p=FALSE)

GOF_results=cbind(round(mf,digits=2),dof,round(pvaluechi,digits=4))


colnames(results)=c('DWMD','SE','z_stat','P_value')

GOF_test=c('Test_value','DOF','P_value','None')
GOF_value=cbind(GOF_results,0)
GOF_results=rbind(GOF_test,GOF_value)
GOF_results

results=rbind(results,GOF_results)
results

namev=c('pi','gamma','tau','theta','lamda','psi','sigma_omega','sigma_csi','sigma_psi','sigma_v','sigma_phi','alpha1')

print(xtable(results, caption="Structural DWMD estimates", digits=3),
      type="latex",file="R:/SharedProjects/Shared2020-070/2016/extend_to_2020/JPE_Replication_log/table7_DWMD.tex",
      caption.placement = "top")
require(foreign)
results_dwmd<-data.frame(cbind(namev,results[1:12,]))
write.dta(results_dwmd,"R:/SharedProjects/Shared2020-070/2016/extend_to_2020/JPE_Replication_dta/R_estimates_DWMD.dta")




